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I. INTRODUCTION 



The underlying idea of the theory of complex sys- 
tems is that complex patterns or structures can be 
generated from simple models; one of the key papers 
leading to this insight was the classic paper of Tur- 
ing pQ on what are now known as reaction-diffusion 
equations. Most of the theoretical studies of pattern 
formation have followed Turing and used partial dif- 
ferential equations (pdes) to specify the model de- 
scribing the system (2). However there is a potential 
problem with this approach: the parameter range for 
which the patterns exist in the model can be very re- 
stricted, in contrast with what is seen in real systems. 
For instance, to observe Turing patterns in simple 
reaction-diffusion systems described by pdes requires 
that the diffusivities of the species be of different or- 
ders [5] [5] . The limited range of parameters for which 
patterns are seen could be attributed to the simplicity 
of the model chosen to describe the process, however 
for systems with an underlying molecular basis an- 
other explanation has recently been put forward [3J- 
[6] . These authors have observed that Turing-like pat- 
terns exist for a much greater range of parameter 
values if the discrete nature of the molecules com- 
prising the system is taken into account. The result- 
ing patterns are sometimes referred to as stochastic 
Turing patterns [5] or quasi- Turing patterns [7J, and 
they may be analyzed using the theory of stochastic 
processes, appropriate given the noise is created as a 
consequence of the discreteness of the system. 

In this paper we investigate a model which not only 
shows Turing and stochastic Turing patterns, but also 
travelling waves and stochastic travelling waves (re- 
ferred to as stochastic waves in the following), the 
latter having the same relation to travelling waves as 
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Turing patterns have to stochastic Turing patterns. 
One interesting aspect of stochastic waves is that they 
can clearly be seen in computer simulations of the 
model. In contrast, while there is unambiguous evi- 
dence for the existence of stochastic Turing patterns 
from, for instance, the form of the power-spectrum 
of the fluctuations, direct visual evidence is less clear 
due to the noisy nature of the patterns. The model 
we study is a Brusselator with a non-local interaction 
term; this choice is largely made on the grounds that 
the model is simple, and so allows the effect to be 
clearly demonstrated. The non-locality seems to be 
an important ingredient in finding stochastic waves, 
similar to the observation that travelling waves dis- 
appear in deterministic models when interactions are 
made more local [5]. Travelling waves have been ob- 
served in chemical reaction systems in [5HTT] , and also 
in other types of population-based systems [T^rUB] . 

The description of patterns through an analy- 
sis of pdes of the reaction-diffusion type is well- 
developed 0, but the analogous stochastic systems 
have received very little attention. The common 
starting point of the work that has been carried 
out [1HS] has been the master equation (continuous- 
time Markov chain) , although the details of the tech- 
niques used to analyze this equation have differed. 
However the whole idea of stochastic patterns, and 
the methods which may be used to analyze them, fol- 
lows on from the idea of stochastic cycles, or quasi- 
cycles, and their detailed study over the last few 
years. Since the essential idea behind stochastic 
patterns can be clarified with reference to the well- 
studied mechanism behind stochastic cycles, we turn 
to their description. 

The context in which stochastic cycles appear is 
in the relation between stochastic individual-based 
models and the corresponding deterministic descrip- 
tions of their dynamics. Interacting many-particle 
systems are typically defined by a set of stochastic 
jTu^Cj^at the microscopic level. Such systems are com- 
inuii In chemistry and in biology [16] . but they are 
also used to model stochastic dynamics in epidemi- 
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ology [T71 [TB], in population dynamics [THl HD] and 
in the social sciences [21) . Frequently in formulating 
descriptions of these systems intrinsic noise, due to 
the discrete nature of the constituents, is treated as 
a negligible perturbation to the dominant determin- 
istic dynamics. It has been known for many years 
that neglecting such fluctuations is not always justi- 
fied, on the contrary intrinsic noise can fundamen- 
tally change the temporal evolution of these systems. 
The activity in this area over the last few years is due 
to the realization that a method imported from statis- 
tical physics — the van Kampen system-size expan- 
sion [16) — can be used to quantitatively understand 
the deviations to the deterministic dynamics caused 
by these stochastic effects. Effectively one uses some 
measure of the inverse system size, for instance its 
inverse volume V~ l , as a perturbation parameter to 
investigate the stochastic dynamics. 

If truncated after the lowest order, the system-size 
expansion provides a systematic path from stochas- 
tic microscopic interacting-particle models to a deter- 
ministic description in terms of differential equations. 
This truncation effectively corresponds to taking the 
thermodynamic limit, i.e. to considering the limit- 
ing case of infinite systems. If carried out to next- 
to-leading order the system-size expansion can suc- 
cessfully characterize stochastic effects in the limit of 
large, but finite systems. Of particular interest is the 
case where the deterministic model has a stable fixed 
point, which is a spiral, so that perturbations away 
from the fixed point decay in an oscillatory manner. 
The effect of the stochasticity is to excite this de- 
caying mode and create sustained oscillations, which 
have an amplitude larger than might naively be ex- 
pected — stochastic amplification 19 . These are the 
stochastic cycles. 

In this description of stochastic cycles no mention 
has been made of the spatial degrees of freedom, and 
indeed most of the studies of stochastic cycles have 
been confined to well-mixed systems where spatial 
effects are ignored. However, in a similar way, intrin- 
sic fluctuations can excite more complicated spatio- 
temporal modes in spatially extended systems, as 
seen for example in a predator-prey model of the 
Volterra type [22]. In this spatially-extended model 
the intrinsic noise leads to spatially uniform tempo- 
ral oscillations, i.e. at zero wavenumber. Stochas- 
tic Turing patterns, just like stochastic cycles, are 
triggered by internal fluctuations, but they occur at 
non-zero wavenumber. To date they have been found 
in a spatial version of the Levin-Segel predator-prey 
model [4 , in the spatial Brusselator model [5] and in 
a model of embryonic patterning [6] . What does not 
change is the ability of the system-size expansion to 
accurately predict the features of the excitations. 

The main purpose of the present paper is to show 
that in addition to uniform spatial oscillations and 
time-independent Turing patterns, intrinsic noise can 



also trigger travelling waves, i.e. phenomena which 
are both spatial and temporal simultaneously. To 
show that such effects may occur we study a variation 
of the Brusselator model with a nonlocal interaction 
term, mediated by an exponential kernel in space. 
We believe that the existence of stochastic travelling 
waves is widespread, and we choose the Brusselator 
model because it is probably the most widely used 
model to illustrate such phenomena, and so is famil- 
iar to a wide number of researchers in the field. A dis- 
advantage is that nonlocal interactions which seem to 
be required to see the effect are less easy to motivate 
in a chemical model such as the Brusselator, but we 
believe that this is outweighed by the advantage of 
using a simple model to illustrate the idea. So while 
several physical mechanisms for nonlocal interactions 
are given in [8], and in biological and social systems 
such effects are far easier to describe and motivate, 
our desire in this paper to avoid specific mechanisms 
and extraneous details, leads us to formulate the non- 
local interaction in a simple and generic way. 

We show that the presence of such a nonlocal term 
promotes travelling waves, and that stochasticity can 
bring about stochastic travelling waves in situations 
where the deterministic system has a stable uni- 
form fixed point. All this is carried out analytically 
through use of the system-size expansion — extended 
to deal with the non-local interaction — and checked 
numerically using the Gillespie algorithm [231 HI] • 

The remainder of the paper is structured as follows. 
In Section [IT] we introduce the model and discuss its 
behavior on the deterministic level. In Section IIIII 
we carry out a detailed analysis of the stochastic dy- 
namics by means of an expansion in the inverse sys- 
tem size. Section |IV| then contains our main results, 
including an analytical characterization of stochas- 
tic waves and confirmation through numerical sim- 
ulations. Our findings are summarized in Section [V] 
where we also discuss possible future lines of research. 
Two appendices contain mathematical details: the 
first a summary of the conditions under which Turing 
and wave instabilities occur and the second aspects 
of the calculations and results from the system-size 
expansion. 



II. MODEL DEFINITION 

The stochastic model is defined as a collection of 
uniform cells, indexed by i, each of which has vol- 
ume V. For convenience these can be thought of as 
cubes in three dimensions, but other regular shapes 
in other dimensions can also be considered. In fact, 
since our main aim here is to illustrate the idea of 
stochastic waves and their characterization, our re- 
sults will largely apply to a one-dimensional system, 
i.e. a chain of cells. Simulations are then less time- 
consuming, but the model still exhibits all features we 
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are interested in studying. Generalization of the ana- 
lytic results to higher-dimensional models is straight- 
forward. Unlike previous work |22j . we will assume 
that the number of cells is infinite, that is, the under- 
lying space is of infinite extent. This avoids technical 
complications when nonlocal interactions are present. 

In every cell molecules of two species, X and 
Y, interact through the reactions of the Brusselator 
model US: 



0, A X h 
Xi — > Yi, 
2X t + Y t 4 3X t 
X,. -A 0„ 



(1) 



The rates at which the reactions occur are denoted 
by a. b, c and d, and X{ and Yi are molecules which 
are in cell i at the time the reaction occurs. The third 
reaction may occur between an A"-molecule in cell i 
and a Y"-molecule in any other cell (with rates de- 
pending on the distance between the two cells), but 
the effect is to reduce the number of Y molecules in 
cell i by 1. This constitutes the nonlocal interac- 
tion; other choices are possible, but this was made on 
grounds of simplicity. The notation used to describe 
the chemical types and the rates in the Brusselator 
model in the literature is not standard, but we follow 
most closely [35] . The number of X and Y molecules 
in cell i will be denoted by 7ij and m, respectively. 
We will also use n and m to represent the spatial 
vectors with components rij and m, respectively. 

The transition rate from the state (n, m) to the 
state (n',m') will be denoted by T(n', m'|n, m), but 
to lighten the notation we will only list the variables 
which have changed in a given reaction. These func- 
tions are found by invoking mass action: 

Ti(m + l,TOj|ni,TOj) = a, 
T 2 {rii - l,mi + l\m,mi) = b y, 



T 3 (m + l,m,i - l\ni,mi) = c A ^ 



j=-oo 



T 4 (rii - l,mi\ni,mi) = d y 



-<r\i-j\ 'J!± 
V ' 



(2) 



where the subscripts on the rates refer to the four 
reactions in Eq. ([lj. These transition rates are as 
in the usual Brusselator model [55], except for the 
third which has the nonlocal character mentioned ear- 
lier. We stress again that it is the influence of the Y 
molecule in cell j that causes the reaction, but that 
the effect is in cell i which is at a distance away. 
The effect is to increase the number of X molecules 
in cell i by one and to decrease the number of Y 
molecules in cell i by one. The form of interaction 
was taken to be exponential because this is a frequent 



choice, but once again on grounds of simplicity. The 
constant a expresses the range of the interaction and 
A is a normalization constant whose choice will be dis- 
cussed later. As a technical aside, the third transition 
rate should also have a factor of 0(mj) — the Heav- 
iside step function — present to prevent the number 
of Y molecules in cell i becoming negative, but given 
we will be looking at a regime far from m, = 0, this 
condition is irrelevant. 

In addition to the reactions given in Eq. 0, "mi- 
gration" reactions which describe molecular diffusion 
from one cell to another, have to be specified. For 
a given cell i, molecules of the two species X and Y 
may diffuse into or diffuse out of a neighboring cell j: 



Xi — > Xi 



X, 



x. 



±>i y Yj . Yj y . 



(3) 



These reactions have transition rates given by 



T 5 (rii - l,rij + l\rii,nj) 
T 6 (ni + l,rij - l\rii,nj) 
Tiirrii — l,m,j + l\mi,mj) 
T 8 (rrii + l,m,j - l\m,i,mj) 



Vz 

Tlj 

H Vz 



(4) 



Here z is the number of nearest neighbors of a cell and 
the index j denotes a nearest neighbor of cell i. For 
the one-dimensional model z = 2 and j E {i— 
It is also worth remarking that we have not imposed a 
fixed limit on the number of molecules permitted in a 
cell (in contrast with some previous work [3 [22] ; this 
is reflected in the fact we use the inverse volume, V^ 1 , 
rather than the inverse total number of molecules, as 
the expansion parameter). 

Since we have assumed that the transition rates 
only depend on the current state of the system, the 
stochastic process is Markov, and so the probability 
distribution for the system being in state (n, m) at 
time t, P n ,m(t), satisfies the master equation [15) : 



dt 



Pn,m(t) = 



T(n,m|n',m')Pn',m'(*) 



T(n',m'|n,m)P n , m (t) 



(5) 



In what follows we will write discrete variables as 
subscripts and continuous variables as arguments of 
functions (such as in P n ,m(i))- The sole exception is 
for the transition rates T(n', m'|n, m), where for the 
sake of clarity we have written the discrete variables 
as arguments of the function T. 

We will discuss the analysis of the master equa- 
tion |5]) using the system-size expansion in the next 
section. For the rest of this section we will obtain 
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the deterministic equation which holds in the limit 
V — > oo directly, without using the system-size ex- 
pansion, and then study the stability of the homoge- 
neous state. 

We begin by defining the concentrations of the X 
and Y molecules in cell i in the infinite volume limit 

by 

Ut)= Hm %^ , Mt) = Jim (6) 

where (. . .) is the mean with respect to the probabil- 
ity distribution P njm (t). Multiplying Eq. ([5| by rii 
and rrii respectively, and using the fact that in the 
deterministic limit the probability distribution is a 
delta- function so that (nf) = (rii) etc, one finds the 
following ordinary differential equations: 



dr 



= a — (b + d)4>i 



ipi-i + a A<j>i 



dr 



ipi = bt 



] = -co 

j = -oo 



A J2 e-'Wfr-i 



+/3 A^i 



(7) 



where r = t/V is a rescaled time and where A/j = 
— 2fi+fi—i is the discrete one-dimensional Lapla- 
cian. 

We choose the normalization constant A so as to 
satisfy 



a J2 e ~ a{jl = L 



(8) 



By doing so, the deterministic equations Q have the 
homogeneous solution: 



a 
d' 



bd 

ac 



(9) 



which are the same as those of the conventional Brus- 
selator model (obtained from our model by replacing 
the non-local interaction by a local term). The ex- 
pression for A can be summed to yield 



A 



1 



(10) 



Eqs. are a set of reaction-diffusion equations of 
the type usually defining the starting point for finding 
Turing patterns. The analysis starts by examining if 
the homogeneous solution (J9j) is unstable to spatially 
inhomogeneous small perturbations. To do this we 
introduce small perturbations 

8<j>i(t) = fait) - <f>* , fyi{t)=il>i(t)-<tl>*, (If) 



into the deterministic equations ^ and keep only 
linear terms in 6<pi(t) and 5ipi(t). This gives 



d 



= — (b + d)Stpi + c(f>* 2 A e-°W Si/ji-j 

j=—oo 

+2 c <f>*il>*8<l>i + a A Sfc, (12) 



dr 



Sfl>i = bS(f>i ~ c<ft* 2 A ^ e-^Sipi-j 

J=—oo 

^c^^Sfa + pASi/ji. (13) 



The structure of these equations makes it clear that 
they will simplify considerably if we go over to a 
Fourier representation. We therefore introduce the 
spatial Fourier transform for the infinite discrete sys- 
tem of cells: 

f(k) = ]T e^ k f s , fj = — dke^ f(k). 

(14) 

Note that A; is a continuous variable which takes val- 
ues in the first Brillouin zone, [0, 2tt\. Fourier trans- 
forming the equations ( 12 ) gives: 



d 

d_ 

8t 



= -(b + d)54> + c(t)* 2 Ae(k)5ij 

+2c(j)*^*S4 l + aAS4>, (15) 



Sip = b 84> - c 4>* 2 A e(fc) 5$ 
-2 c 4>* iji* 54> + ft A 6$, 



(16) 



where S(f> — S(f>(k,T) and similarly for Sip. The two 
functions e(fc) and A = A(fc) are respectively the 
Fourier transform of the exponential function and of 
the Laplacian, i.e. 

oo 

Af(k)= J2 (fi+i- 2/,: + /i-i) = 

j=-oo 

= 2 [cos(fc) - 1] f(k) => A = 2 [cos(fc) - 1] , 



e(fc)= £ e 



sinh(er) 



j=-oo 



cosh(cr) — cos(fc) 



(17) 



The system (151 may be written in a more compact 
form: 



(18) 



where 



J*(k) = 

-(b + d) + 2c0*V* +« A cAf 2 e(fc) \ 
b-2c(p*ip* -c A cj)* 2 e(fc) + p A J ' 

(19) 
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The eigenvalues of the Jacobian matrix (19), Ai(fc) 
and \2{k), yield information about whether perturb- 
ing the homogeneous solution leads to pattern forma- 
tion. If both Ai and A2 have negative real part (i.e. 
Re[Aj(fc)] < 0, V7c, i — 1,2), the homogeneous state 
is stable: every perturbation will eventually die out 
and no pattern will develop. If, on the other hand, 
there is an eigenvalue at a non-zero k with a positive 
real part, then a spatially modulated instability oc- 
curs: a perturbation will grow in magnitude taking 
the system from the homogeneous state to one with 
wavenumber defined by k. This growth will eventu- 
ally be saturated by the non-linear terms leading to 
a pattern of characteristic wave number k. 

This linear analysis of the homogeneous state is 
also able to determine whether the resulting pattern 
is steady or oscillatory, by looking at the imaginary 
part of the eigenvalues u>i = Im[A,]. Steady patterns 
correspond to Im[A-j(fc)] = 0, for all unstable modes k, 
the case in which the instability is called a "Turing 
instability". When Im[A,(fc)] 7^ 0, for an unstable 
mode at a non-zero k, the system is said to undergo a 
"wave instability" as the resulting pattern will consist 
of travelling waves [2]. 

When multiple instabilities occur simultaneously, 
say of the Turing type for some k\ and the wave type 
for some k 2 , it becomes harder to predict if the final 
pattern will be steady or oscillatory. However, for 
model parameters sufficiently close to the region in 
which the homogeneous state is stable, only a single 
instability occurs. We will therefore study the insta- 
bilities at the border of the stable region, saying that 
there is a Turing instability (resp. wave instability) 
when the instability near the border is of the Tur- 
ing type (resp. wave type). The criteria we use to 
determine the stability boundaries are discussed in 
Appendix A. 



The above analysis when applied to Eq. ( 19 1 is dis 



cussed in Section |IV| However, as explained in the 
Introduction we wish to go beyond this deterministic 
approximation, and examine patterns which emerge 
from stochastic effects which are a consequence of the 
underlying discreteness of the system. We therefore 
now go on to discuss the stochastic analysis of the 
model, before collecting results on the determinstic 



For example, the term in the master equation which 
involves the first reaction in Eq. pj , and so corre- 
sponds to the function Tj in Eq. (J2T would be given 
by 

i 

= E( e ^- 1 ) aP n,m(*). (21) 



The master equation rewritten in this way, with all 
eight terms present, is given by Eq. (Bl ) in Appendix 
B. 

We will now carry out the system-size expansion. It 
is important to realize that, for our model, this is an 
expansion in powers of V -1 / 2 , where V is the volume 
per cell. The expansion therefore captures stochastic 
effects at large, but finite cell volumes. This limit 
should not be confused with the limit of a infinite 
number of cells. When we use the term 'system-size 
expansion', we always refer to an expansion in the size 
(volume) of a single cell. Even if the volume per cell 
is finite, the total number of molecules in the system 
(summed over all cells) may well be infinite, if the 
system is composed of an infinite number of cells. 

The system-size expansion itself involves making 
the time-dependent change of variables (n, m) n- 



n ^ V4>{t) + Wt, m i-> Vif){t) + Wt), (22) 

where <p(t) and i/)(t) are two time-dependent vectors 
defined by Eq. each of these vectors having as 
many components as there are different cells in the 
system. From the leading-order term in the expan- 
sion one finds that these quantities satisfy the deter- 
ministic equation ([7]). 

We now change the degrees of freedom of the 
stochastic system to £ and rj and consider the prob- 
ability distribution in terms of these variables as 



n(£,T7,t) = p n>m (t). 



(23) 



and stochastic regimes in Sec. IV 



The left-hand side of the master equation ^ has the 
form 116 



III. STOCHASTIC ANALYSIS 

The analysis of the model, without making the de- 
terministic approximation, begins by writing down 
the master equation ^ in a form which is more 
amenable to application of the system-size expansion. 
This is done by introducing step-operators |16j : 



ef ,i/(n, m) =/(..., n t ± 1, . . . , m), 
4,i/( n > m ) = /( n ; • • • ,Wi ± 1, . . .)• 



(20) 



dP 
~dt 



d t H - W V € n • d t cf> - Vv \/ v U ■ d t ip, (24) 



where d t = d/dt. 

To determine the nature of the expansion of the 
right-hand side of the master equation we begin with 
the form (Bl) given in terms of the step operators. 



This is because the step ope rato rs have a natural ex- 
pansion in V - s given in Eq. (B2 ) of Appendix B. The 
introduction of a rescaled time, r = t/V brings the 
master equation into the general form (see Appendix 
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<9 T II - -= (V € n • d T cf> - Vr,U ■ d T 1p) 



V 

[ y [f(<p, ik) ■ v € + g(</», t/>) • v„] + ^ n. 



B) 
1 

V 



(25) 

where L is a linear operator containing various deriva- 
tives in 7/ and £ and / and g are functions of (p and 
i>. 

It is now possible to match terms on both sides 
of the transformed master equation (25 1. The order 
1/y/V contributions lead to 



d_ 
d_ 



<pi = Mcp,ip), 

ipi = gi(cp,ip), 



(26) 



which are just the deterministic equations for the con- 
centrations in cell i derived in a more direct fashion 
earlier, and given explicitly by Eq. (JT]). Matching the 
order 1/V contributions leads to an equation for the 
probability distribution II which describes the fluctu- 
ations: 



P(^r,,t) = Ln(ti?,t). 



(27) 



We have analyzed the leading order result (261 in 
Section [TTJ our aim in this Section is to study the 
stochastic corrections given by Eq. ( 27 ) . The explicit 



form of the operator L given in Appendix B shows 
Eq. (27) to be a Fokker-Planck equation describing a 
linear stochastic process. It is more convenient for our 
purposes to use the equivalent description of the pro- 
cess in terms of a linear stochastic differential equa- 
tion, or Langevin equation, since this allows us to 
take the spatial and temporal Fourier transforms. If 
we carry out a Fourier transform only with respect to 
the spatial variable, one has the form [57J HE] 



dr 



i(C) + £(fc,r) 



(28) 



where fi{k, r) is a Gaussian noise with zero mean and 
correlator 

(A(fc,r)/i T (A:',T')) =B{k)2-K5{k-k')8{T-T J ), (29) 

and where we have introduced the convenient nota- 
tion C = (£, »7) • 

The explicit forms for A and B are given in Ap- 
pendix B. Since the stochastic process is linear, A 
is a linear function of £ and B is independent of it. 
Both A and B are explicit functions of r through 
their dependence on the solutions of the deterministic 
equations 4>{t) and ip(r). However we are interested 
in fluctuations about the homogeneous state, and so 
we may take its values for these solutions. In this 



case both A* and B* lose their explicit time depen- 
dence. If the noise term p, is omitted from Eq. ( 28 ) 
(effectively by taking B — 0; see Eq. (29)), then £ is 



nothing else than a deterministic perturbation of the 
deterministic dynamics. So it is not surprising that 
A* is related to the elements of the Jacobian in the 
homogeneous state given by Eq. (19): 



A*Jk, 



J2Jr*MUk,r) 



(30) 



s = l 



Of course A retains an implicit dependence on time 
through £(k,r). The form for B{k) is given by 
Eq. (JbTT ) . The Langevin equation (28) can now 



be solved by taking the temporal Fourier transform 
and determining the Fourier-transformed fluctuations 
£(k,u) and fj(k,u) [22] . 

For parameter values for which the homogeneous 
state is stable no pattern arises in the deterministic 
description. However, as discussed in the Introduc- 
tion, we can ask if in this region of parameters the 
spatial system exhibits ordered structures once fluc- 
tuations are taken into account. To probe this possi- 
ble fluctuationally-induced order, a useful tool is the 
power spectrum of the fluctuations about the homo- 
geneous state, defined by 



Px(k,w) 



Py(k,u) = (\fj{k,w)\ 



(31) 

In absence of order the spectra will show an almost 
flat profile. If instead some type of order is present, 
the power spectra will typically have a characteristic 
peak. The position of the peak in combined (fc, uo)- 
space determines the type of structure that is present. 
For example, a global oscillation in time will corre- 
spond to a peak of P{k, uj) at a non-zero value of uj 
and at k = 0, whereas the power spectrum peaks at 
uj = and at a non-zero value of k for stochastic Tur- 
ing patterns. As we are looking for stochastic waves 
we shall seek parameter values for which the power 
spectra display a peak at values of (fc, uj) where both 
k 7^ and uj ^ 0. The spectra are found both through 
simulations and analytically from the form of £(fc, uj) 
and ij(k,uj), derived from Eq. (28). The analytical 



calculations required to determine the power spectra 
are similar to those discussed in [25], but with dif- 
ferent forms for the matrices J*(k) and B*{k). One 
finds 



P x (oj,k) 
Py(u,k) 



Cx + B\ uj 2 

{uj 2 - n 2 ) + r 2 uj 2 

C Y + B\ 2 u 2 
JuV 2 



n 2 ) 



r 2 , 



(32) 



where fi = ^det J*(k), T = -tr J*(k) and where 
Cx(k) and Cy(k) are given in Appendix B. Expres- 
sions ( 32 ) constitute a full analytical characterization 
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Figure 1 : Phase diagram in the (c, /3) plane for the de- 
terministic equations Q obtained with a linear analysis 
of the homogeneous state (|9| for a = d = a = X, <r = 2 
and 6 = 10. The model exhibits a phase in which the ho- 
mogeneous fixed point is stable, along with phases with 
Turing patterns and travelling-waves. The two lines in the 
diagram, j3w(c) and /3t(c), indicate the onset of these in- 
stabilities. Far from these lines simultaneous instabilities 
may occur (e.g. travelling waves and Turing patterns). 



of the spatio-temporal power spectra of fluctuations, 
or equivalently of their spatio-temporal correlation 
properties. This completes the mathematical analy- 
sis within the van Kampen formalism. The form of 
the spectra obtained, and their interpretation will be 
discussed in the next section. 



IV. RESULTS 

The main purpose of this paper is to investigate 
how intrinsic fluctuations can bring about stochastic 
waves. Before we turn to the stochastic dynamics it 
is convenient to characterize the behavior of the non- 
local Brusselator model in the deterministic approxi- 
mation. Having derived analytical expressions for the 
homogeneous state of the deterministic system (see 
Eq. (|9|) , as well as for the corresponding Jacobian, 
Eq. |T9j ), it is straightforward to obtain the stability 
properties of the system as a function of the model 
parameters. This is discussed at the end of Sec. [TT] 
and in Appendix A. Given the comparatively large 
number of parameters, we focus on a representative 
selection of phase diagrams. 

The properties of the deterministic dynamics at 
varying values of the parameters /3 and c (and keeping 
all other parameters fixed) are illustrated in Fig. [I] 
For a fixed value of c we find a phase at intermedi- 
ate values of j3 € [ftw , Pt\ in which the homogeneous 
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state 



waves 







1 



cr 



Figure 2: Phase diagram for the deterministic dynamics, 
obtained from Eqs. (7| for a = d = a = 1, P — 0.1 
and b — 10. The solid line marks the onset of a wave 
instability, c = cw (cr) . The diagram illustrates the role of 
the nonlocal interaction in generating a wave instability. 



state is stable against fluctuations of any wave num- 
ber. At a critical value /3 = /3t(c) a Turing instabil- 
ity sets in, i.e. an unstable mode occurs at a nonzero 
wave number, with the corresponding eigenvalue be- 
ing real. At values of j3 lower than some second crit- 
ical value, /3w(c), the instability occurs again at a 
non-zero wave number, but now the corresponding 
eigenvalue is complex, indicating a wave instability. 
We have thus established that the nonlocal Brusse- 
lator model exhibits Turing instabilities, as well as 
travelling- wave instabilities in the deterministic limit. 

In order to illustrate the role of the nonlocal inter- 
action term we show a second phase diagram, now 
in the (cr, c)-plane, in Fig. [2] Recall that the model 
parameter cr characterizes the range of the nonlocal 
interaction: for small values of cr the interaction ker- 
nel in Eqs. ^ decays slowly with distance, and the 
interaction is therefore long range. For large a the 
interaction range is small, in the limit a — > oo one 
recovers the standard Brusselator model with purely 
local interactions. This is clear from Eq. |2]); the only 
term in the sum which is independent of a is the j = 
term, all the others exponentially decay with cr, and 
so vanish as a — > oo. Thus the sum tends to rrii/V 
as a — > oo. 

As seen in Fig. [2j the long-range nature of the in- 
teractions can promote the occurrence of travelling 
waves in the deterministic system. Varying a at fixed 
value of c (and all other model parameters fixed as 
indicated in the caption of Fig. [2]) one finds that the 
homogeneous state is stable at large values of cr, i.e. 
when interactions are localized. As a is reduced (the 
interaction range is increased) a wave instability sets 
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Figure 3: (Color online) Spatio-temporal dynamics of 
the concentration of Y-molecules in the stochastic waves 
regime. The number of particles, obtained from a sin- 
gle run of the stochastic dynamics, is reported on a scale 
of shades (with lighter shades indicating a larger number 
of Y-particles per cell, and dark shades representing low 
concentrations). Stochastic waves are seen as diagonal 
structures in the figure, representing both right-moving 
and left-moving waves. The system considered is a one- 
dimensional ring of 20 cells, each of which has a volume 
V = 5000. Parameter values are a = d = a = 1, <J = 2, 
c= b = 10 and f3 = 0.1. 



We now turn to the stochastic system and show the 
spatio-temporal behavior of the concentration of Y- 
molecules in Fig. [3] It is important to stress that we 
have chosen values of the model parameters such that 
the deterministic system has a stable homogeneous 
state. More specifically, for the parameters chosen in 
Fig. [3j the system is in the stable phase of Fig. [2j but 
near the line along which a wave instability occurs 
in the deterministic model. As demonstrated by the 
diagonal structures in the space-time representation 
of Fig. [3j the stochastic system displays travelling 
stochastic waves for this choice of parameters, even 
if such waves are absent in the deterministic system. 
In order to illustrate the mechanism producing these 
stochastic waves we plot the dispersion relations of 
modes near the wave instability of the deterministic 
system in Fig. [4] (model parameters other than c are 
as in the simulations shown in Fig. [3]). For c > 9.35 
all modes are stable, and the least stable mode (i.e. 
the mode for which the real part of the corresponding 
eigenvalue is maximal among all modes) has an eigen- 
value with non-zero imaginary part and negative real 
part. At c w 9.35 the wave instability line is crossed, 
one mode is now marginally stable, and the corre- 
sponding eigenvalue is purely imaginary. Crucially, 




Figure 4: Real part (solid line) and imaginary part (inset, 
dashed line) of A(fc), one of eigenvalues of the Jacobian 
matrix (19 1 (the other is its complex conjugate). Model 
parameters are a = d = a = 1, a = 2, b — 10 and (3 = 0.1. 
Using symmetry, we restrict the range of k to a half of the 
Brillouin zone, [0,7r]. The onset of the travelling wave 
instability occurs at c ~ 9.35. 
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Figure 5: Upper panel: Temporal evolution of the Y- 
concentration in selected cell in the stochastic waves 
regime. Data is obtained from one run of the stochas- 
tic dynamics, with parameter values as in Fig. [3] Lower 
panel: Snapshot of the concentration of y-molecules as a 
function of position at a fixed time t. 



the instability occurs at a non-zero wave number, and 
with a non-zero imaginary part of the corresponding 
eigenvalue (see inset). For c < 9.35 the deterministic 
system has unstable modes and it exhibits travelling 
waves. 

The stochastic simulations of Fig. [3] are carried out 
at c = 10. Here all modes are stable in the deter- 
ministic system, the least stable mode has a non-zero 
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wavenumber and a complex eigenvalue. In the ab- 
sence of intrinsic stochasticity the system would con- 
verge to the deterministic homogeneous state. Fluc- 
tuations, due to demographic noise at finite cell 
volumes, however constantly cause random motion 
about the homogeneous state. At large, but finite 
cell volumes, a linear approximation is in order and 
the fluctuations can be decomposed into their Fourier 
modes. The fluctuations corresponding to the least 
stable mode decay the slowest, on a time scale set by 
the real part of the corresponding eigenvalue. This 
effect, along with fluctuations persistently occurring, 
conspire due to the mechanism now known as co- 
herent stochastic amplification [19 . Modes which 
are stable in the deterministic model can effectively 
be excited by random fluctuations, and are observed 
with sizeable amplitude in the stochastic system. In 
our model this amplitude scales as V~ 1//2 in the cell 
volume, but the prefactor multiplying this factor can 
be significant, resulting in stochastic waves with ap- 
preciable amplitude even at large cell volumes. 

To illustrate the occurrence of stochastic waves fur- 
ther, we show a time-series of the (re-scaled) con- 
centration of Y molecules in a fixed cell in the up- 
per panel of Fig. [5] Coherent stochastic oscillations 
are clearly visible, their amplitude is seen to scale 
as V -1 / 2 with the cell volume. Plotting the time 
series of the Y-concentration in one cell we have ef- 
fectively disregarded the spatial structure of the sys- 
tem, and focused instead on the temporal dynamics 
only. In this sense the upper panel of Fig. [5] is simi- 
lar to observations of amplified stochastic oscillations 
in the time series of non-spatial systems (see for ex- 
ample [T9]). The lower panel of Fig. [5] shows the 
concentration of Y molecules as a function of posi- 
tion at a fixed time. Spatial modulations, again with 
an amplitude scaling as V~ 1 / 2 , are seen, even though 
of a lower coherence than the temporal oscillations 
shown in the upper panel. In this lower panel we 
have effectively disregarded the temporal dynamics 
of the system, and have instead focused on its spatial 
character. The modulations of the concentration of 
Y molecules as a function of position therefore consti- 
tutes the analog of stochastically-driven Turing pat- 
terns reported in [3H5] . The novelty of our model is 
that it combines the spatial and the temporal aspect; 
the stochastic waves seen in the nonlocal Brusselator 
model are noise-driven patterns with structure both 
as a function of position and as a function of time. 



The theory developed in Sec. Ill allows us to char- 
acterize the observed stochastic waves further. Based 



on the results of Eqs. (32) we are able to predict the 



power spectra of spatio-temporal fluctuations about 
the homogeneous state for any choice of model pa- 
rameters. Stochastic waves are to be expected when- 
ever the highest peak of these power spectra is found 
simultaneously at a non-zero wavenumber k (indicat- 
ing non-trivial spatial structure) and at a non-zero 





Figure 6: (Color online) Power spectrum of the fluctua- 
tions of th e Y species obtained analytically (upper panel) 
from Eq. ( 32 1 and numerically (lower panel) by simulat- 



ing the stochastic process using the Gillespie algorithm. 
The agreement between the two power spectra is clearly 
very good. The system is in the stochastic wave regime 
with a = d = a = l, o = 2, c = 6= 10 and /3 = 0.05. The 
spectra show a peaked profile which corresponds to spatio- 
temporal organization despite the deterministic predic- 
tion of a stable homogeneous state. The numerical spec- 
trum is obtained by averaging 200 realizations of a finite 
system of 28 cells each of which has volume V = 1500. 



angular frequency uj (indicating a complex eigenvalue, 
and hence oscillatory behavior as a function of time). 
An example of a power spectrum with these proper- 
ties, obtained from the theory of Sec. 



Ill is shown 



in the upper panel of Fig. [6j for comparison we show 
measurements from direct numerical simulations in 
the lower panel. As seen from this Figure, the agree- 
ment between theory and simulations is excellent. 



V. CONCLUSION 

This paper has been concerned with the investiga- 
tion of stochastic effects in a nonlocal extension of 
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the Brusselator model. An analysis of this model on 
the deterministic level reveals that nonlocal interac- 
tions can promote the occurrence of travelling wave 
instabilities, similar to what has been seen before in 
other chemical reaction models with long-range in- 
teractions [8|. As the main result of our work we 
show that the nonlocal Brusselator model can also 
exhibit travelling waves driven by internal fluctua- 
tions in parameter regimes in which the deterministic 
system converges to a homogeneous state. Based on 
a stochastic formulation of the model in terms of a 
chemical master equation, and a subsequent expan- 
sion in the inverse system size, we derived analyti- 
cal expressions for the power spectra of these spatio- 
temporal patterns, in excellent agreement with direct 
numerical simulations. 

These findings extend previous results on noise- 
driven instabilities. In |19| non-spatial systems were 
considered, and it was shown that intrinsic fluctua- 
tions can generate coherent stochastic oscillations for 
parameters for which the deterministic system spi- 
rals into a fixed point. The work of j4j-[6] instead 
focused on spatial systems, and it was shown that in- 
trinsic noise can generate stochastic Turing patterns, 
i.e. spatial structures with a constant amplitude in 
time. Our model combines both aspects, and pro- 
duces stochastic patterns with a full spatio-temporal 
dynamics. 

These stochastic waves are seen in the power spec- 
tra of fluctuations, computed from the theoretical 
approach, as isolated peaks at non-zero wave num- 
ber and non-zero angular frequency. While in the 
past it has been rather difficult to observe spatial 
stochastic patterns directly, and while most of the 
previous work was limited to an indirect identifi- 
cation in Fourier space, we have also been able to 
obtain direct visual confirmation of the stochastic 
waves in the Brusselator model with nonlocal interac- 
tion. Criteria which distinguish stochastic cycles or 
stochastic patterns from their deterministic analogs 
(limit cycles and Turing-like patterns) have been pro- 
posed |H 13 [55]. However these are not applicable 
to systems with nonlocal interactions of the type we 
have considered here, due to the extra fc-dependence 
which comes about because of the nonlocality. It 
would be interesting to devise criteria which encom- 
pass nonlocal models as well. 



The purpose of the current work is to illustrate a 
generic phenomenon which is expected to occur in 
a wide class of systems; the Brusselator model was 
chosen to illustrate the basic idea because it is sim- 
ple and widely studied. The effects of noise on trav- 
elling waves has been investigated in the past, but 
most of these studies depended on the properties of 
the specific system under consideration. For exam- 
ple, a simple reaction diffusion equation displaying 
travelling waves is the Fisher equation [30]. For this 
system it is known that noise may alter the prop- 
erties of the waves (see for instance [31]) and may 
also be responsible of the emergence of noise-driven 
travelling waves |32j . However, these waves are dif- 
ferent to those we have discussed here, as they do 
not arise from an underlying unstable homogeneous 
state. A related point is that the Fisher equation has 
only one species, whereas we need at least two inter- 
acting species, since we require complex eigenvalues 
to trigger the wave instability [5] . 

We expect that the mechanism of coherent am- 
plification, applies to more complicated linear insta- 
bilities as well [3] and that the concept of stochas- 
tic waves is relevant in other models with travel- 
ling fronts. We also expect that the analytical for- 
malism we have developed here for spatial systems 
with nonlocal interaction can successfully be em- 
ployed to study broader classes of individual-based 
models. This may include models of the spread of 
epidemics, in which infection can occur at a distance; 
deterministic models of such processes have been con- 
sidered [331 [31] ■ As stochastic patterns are found in 
more and more model systems, we expect the search 
for them in real systems to intensify and the un- 
derstanding of the underlying causes of pattern for- 
mation in physical, biological and social systems to 
broaden. 
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Appendix A: Turing and Wave Instabilities 

In this Appendix we derive the conditions which 
allow the onset of Turing or wave instabilities to 
be located. Such criteria exist in the literature [T]- 
[3], but most of them rely on the standard reaction- 
diffusion paradigm and will not be applicable in our 
case which includes a nonlocal kernel, and so an extra 
fc-dependence over and above that coming from the 
diffusion terms. We need therefore a more general 
approach, described below, which closely follows |35j . 

We start by defining the region of parameter space 
in which the homogeneous state is stable. Its borders 
delimitate the instabilities whose type can be deter- 
mined from the analysis below. The stability con- 
dition, Re[Ai(fc)] < for all fc, can be conveniently 
rewritten using the trace and determinant if the Ja- 
cobian is a 2 x 2 matrix as 

dctjT(fc)>0, tr.7*(fc)<0, Vfc. (Al) 

The stability region is the set of parameters 
in which the above inequalities hold. Plotting 
det J* (fc c ) against tr J* (k c ) we see that we may leave 
the stability region by violating one of these inequal- 
ities. That is, when 

• There exists a k c ^ such as det J*(kc) = 
whereas trj*(k) < Vfc. 

• There exists a kc ^ such as tr J* (kc) = 
whereas det J* (k) > Vfc. 

It is also possible that determinant and trace be- 
come simultaneously zero, but this is a degenerate 
case which we do not consider here. 

Now recall that the eigenvalues of the Jacobian are 
given by 

Ai, 2 = \ ftr J* ± (tr J*) 2 ~ 4 det J*^j , (A2) 

from which the imaginary part of the eigenvalues may 
be found. From the discussion at the end of Sec.Hllwe 
can see that the above conditions defining the bound- 
aries of the stability region correspond respectively to 
a Turing instability (the former) and to a wave insta- 
bility (the latter). 

The stability conditions as given are not so conve- 
nient to deal with directly, because of the presence 
of inequalities which must be solved for every fc. To 
overcome this, we suppose that trj*(k) has a global 
maximum at ku and det J*(k) a global minimum at 
fc m , a hypothesis that will be discussed further below. 
In this case the two conditions may be rewritten as: 

• (Turing instability) det J* (k m ) = and 
tiJ*(k M ) < 0. 



• (Wave instability) trJ*(kM) 
detj*(k m ) > 0. 



= and 
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These are the conditions we have used to obtain Fig.[T] 
and Fig. [2j 

We can now check that the specific forms of 
tvj*(k) and det J*{k) in our model have the re- 
quired extrema. Finding the extremal points of 
trJ'*(fcAf) can easily be achieved analytically; for 
det J* [k m ) it is a little more difficult. However check- 
ing the existence of a global maximum or minimum 
numerically is straightforward, and we have verified 
that for the range of parameters of interest to us in 
this paper such extrema always exist and moreover 
give the boundaries shown in Fig. [T] and Fig. [2j 



Appendix B: The van Kampen System-Size 
Expansion 

A description of the general structure and method- 
ology behind the system-size expansion, as applied to 
the system under consideration, is given in Section 
III In this Appendix we give some of the technical 
details which would otherwise disrupt the flow of the 
arguments in the main text. 

The application of the method is facilitated by 
writing down the master equation ^ in terms of the 
step operators (20 1. An example is given in Eq. (21 1 



for the first reaction of the set of reactions given by 
Eq. ([!]). When all eight reactions are included the 
master equation is given by 

d „ 

oo 

= E ( e X,i~ 1 ) T ( n i + 1 ' m i\ n i>'ini) 
i— — oo _ 

+( e x,t ~ !) T ( n i — l,mi, \rii, to^ 



E ( :: 'a. 

j'e{i-l,i+l} 



sions [IT)] : 



-X,i 



1± 



V 



-YA 



1± 



V 



(B2) 



This together with the substitution of the ansatz ( |22[ ) 
into the transition rates ([2]), and the replacement of 
Pn,m(t) by IL(£,T),t), allows us to expand the right- 
hand side in powers of V" 3 . 

Equating the left- and right-hand sides of the mas- 
ter equation gives the general form 

= [f{<t>, i>) ■ v c + g{4>, V) • v„] + y\ n, 



(B3) 

where L is a linear operator containing various deriva- 
tives in rj and £ and / and g are functions of and 
After the introduction of a re-scaled time r = t/V, 
Eq. (251 of the main text is obtained. This is now in 



a form where the various terms on both sides of the 
equation can be balanced. 

Equating the leading order terms in Eq. ( [25] ) gives 
equations whose general structure is displayed in 
Eq. (26) and whose specific form in Eq. (|7). Equat- 
ing the next-to-leading order terms gives an equation 
with general structure (27) and specific form 



1) T{m + 1, mi - l\m, rrn) 

1) T(n,; — l,rij + l\rii,nj) 



i— — oo \ r— 1 

+\ E E d <sArA B rs,ijn]) ,(B4) 

r,s— 1 j=i — 1 / 

where for convenience we have introduced the nota- 
tion £i = £ and (2 = V- The precise forms of A r ,i and 
of B rs Aj will be discussed below, but as mentioned 
in the main text it is easier for our purposes to work 
with the equivalent Langevin equation [271 128] 



X,j "X,i 



1) T{m + l,n 3 - l\m,rij) 



dr 



A-(C) + Mi(r), 



(B5) 



+( e ?j e Y,i _ !) T( m t + 1- m j - M m n m j) 
+( e Y,i e Y.j ~ !) T ( m i - 1) m i + M m i, mj) 

XP„,m(*). 



where is a Gaussian noise with zero mean and 
correlator 



(m(T)tf(7>))=B ij S{T-T J ). 



(B6) 



(Bl) 

Wc arc then able to take the Fourier transform of this 



The fundamental ansatz of the system-size expan- 
sion is given by Eq. (22) in the main text. It leads 
to the expression (24) for the left-hand side of the 



equation in both space and time to give Eq. ( 28 ) 



The Langevin equation ( 28 1 is defined by two ma- 
tricies: the drift matrix and the diffusion matrix. 
Since the drift matrix is identical to the Jacobian of 



master equation. The right-hand side can be evalu- 
ated by observing that the step operators (20) can 



the system, as presented in Eq. ( 30 ), we will only give 



be expanded in the inverse of the square root of the 
system size, V~ 2 , giving rise to the following expres- 



the expression for the diffusion matrix. 
The diffusion matrix B has elements B r 



where 



r and s index the species and i and j the cell. In 
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the following the expressions are given for each r and 
s and for a given cell i. The only non-zero values 
of B rSt ij occur when j = i — 1, i or i + 1, and these 
are given respectively as the first, second and third 
entries of a row vector: 

= f-a^i + a+(i> + d)0i 

OO 

+cA ^ e- <r l'l#Vi-i+a(&-i + 2& + &+i), 

J=-0O 



^+1 



= 10, -I 



B*2A = 



OO 

j=-oo 

(B7) 



Evaluating them in the homogeneous state gives 
R* — 



2aa ^ ^ab 4aa 
d ' d d ' 



2aa 
~d~ 



n* — K* — 
"12, i — °1\A — 

Bh,i -- 



0, -« 



2M/3 2a& 46d/3 26d/3 
ac ' d ac ' etc 



(B8) 



The structure of <B* Si jj can be seen from Eq. (B8) 
to be 



B?.,ii = $ ) *W,o + && ) Vil,i. ( B9 ) 



where the two matrices and can be read off 



from Eq. ( B8 1 . It is then straightforward [22] to cal- 
culate the spatial Fourier transform, B rs = B rs (k), of 
the matrices (B8) with respect to the variable i — j: 



B r * s (fc) = (&(° s )+2&«)+&WA, 



(BIO) 



where A is given by Eq. (17 1. The explicit forms are 



B\ 2 (k) - B* 21 (k) 



2ab 



2ab 2bdf3 - 



(Bll) 



The power spectra of the stochastic oscillations, de- 



fined by Eq. (31), have the form (132]) . The functions 



Cx and Cy are defined by 



Cx{k) 



SiiW J22 (k) 2 - 2 B* 12 (k) J* 2 (k) J* 2 (k) 
+B* 22 (k)J 1 * 2 (k) 2 , 



C Y (k) = B* 22 {k)J* 1 {k) 2 -2B* 12 {k)J 2 \{k)J* 1 {k) 



+J3* n (k) J 2 \{kf 



(B12) 



